Diffusion in sparse networks: linear to semi-linear crossover 



Yaron de Leeuw and Doron Cohen 
Department of Physics, Ben Gurion University of the Negev, Beer Sheva 84105, Israel 

We consider random networks whose dynamics is described by a rate equation, with transition 
rates w„m that form a symmetric matrix. The long time evolution of the system is characterized by 
a diffusion coefficient D. In one dimension it is well known that D can display an abrupt percolation- 
like transition from diffusion (D > 0) to sub-diffusion (D — 0). A question arises whether such a 
transition happens in higher dimensions. Numerically D can be evaluated using a resistor network 
calculation, or optionally it can be deduced from the spectral properties of the system. Contrary 
to a recent expectation that is based on a renormalization-group analysis, we deduce that D is 
finite; suggest an "effective-range-hopping" procedure to evaluate it; and contrast the results with 
the linear estimate. The same approach is useful for the analysis of networks that are described by 
quasi-one-dimensional sparse banded matrices. 



I. INTRODUCTION 

The study of network systems is of interest in diverse 
fields of Mathematics, Physics, and Computer and Life 
Sciences. Commonly a network is described by a sym- 
metric matrix that consists of real non-negative elements, 
e.g. the adjacency matrix, or the link probability ma- 
trix, that have unique spectral properties [H [2] . Physi- 
cally motivated, in this work we consider d-dimensional 
network systems, whose dynamics is described by a rate 
equation 



dPn 
dt 



^ ^ tVn rnP m 



(1) 



The off-diagonal elements of w are the transition rates, 
while the diagonal elements are the decay rates 
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We assume a symmetric matrix and write schematically 



w = matrix{w„ m } 



(3) 



In some sense, one can regard to as a discrete Laplacian 
that is associated with the network. Clearly the physi- 
cal problem is related to the study of random walk in a 
disordered environment [3HS] ■ 

For presentation purposes we regard the nodes of the 
network as sites, each having a location x n . By con- 
struction, we assume that the transition rates w nm are 
given by the expression WQe~ enm B(x n — x m ), where B(r) 
describes the systematic dependence of the coupling on 
the distance between the sites, and e is a random variable 
that might represent, say, the activation energy that is 
required to make a transition. Consequently the network 
is characterized by two functions: 

w(r, e) = w e~ e B(r) (4) 
p(r, e) = local density of sites (5) 

The latter is defined as the density of sites in (r, e) space, 
relative to some initial site. Obviously the functional 



dependence of this density on r is affected by the dimen- 
sionality of the network. 

Sparsity.- Our interest is focused on "sparse" net- 
works. This means that the transition rates between 
neighboring sites are log-wide distributed as in "glassy" 
systems. These rates span several orders of magnitudes 
as determined by the dispersion of r or by the dispersion 
of e. In particular (but not exclusively) we are interested 
in a random site model where the rates depend expo- 
nentially on the distance between randomly distributed 
sites, namely B(r) = exp(— r/£). In this particular case 
one can characterize the sparsity by the parameter 



C/ro 



(6) 



where tq is the average distance between neighboring 
sites. We refer to such networks as "sparse" if s <C 1. 

Sparsity vs percolation. - The problem that we consider 
is a variant of the percolation problem [6]: Instead of 
considering a bi- modal distribution ( "zeros" and "ones" ) 
we consider a log- wide distribution of rates [7] , for which 
the median is much smaller than the mean value. We call 
such a network "sparse" (with quotation marks) because 
the large elements constitute a minority. 

Sparsity vs disorder.- While the standard "percola- 
tion" problem can be regarded as the outcome of extreme 
"sparsity" , the latter can be regarded as arising from an 
extreme "disorder" . Accordingly, the model that we are 
considering is a close relative of the Anderson localization 
problem, and therefore we shall dedicate some discussion 
to clarify the relation. 

Physical context.- The model that we address is re- 
lated and motivated by various physical problems, for 
example: phonon propagation in disordered solids [SJ- 
H"0] ; Mott hopping conductance [Tllfl5] ; transport in 
oil reservoirs [HI [T7] ; conductance of ballistic rings [TB] ; 
and energy absorption by trapped atoms |19j . Optionally 
these models can be fabricated by combining oscillators: 
say mechanical springs or electrical resistor-capacitor el- 
ements. In all these examples the issue is to understand 
how the transport is affected by the sparsity of a network. 
If the rates are induced by a driving source, this issue can 
be phrased as going beyond the familiar framework of lin- 
ear response theory (LRT), as explained below. 
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Diffusion and subdiffusion.- Our interest is focused 
on the diffusion coefficient D that characterizes the long 
time dynamics of a spreading distribution. The simplest 
way to define it, as in standard textbooks, is via the 
variance S(t) = (r 2 ).. Namely, 



D = {2d)-' lim 



S(t) 



t 



(7) 



Optionally it can be defined or deduced from the decay 
of the survival probability V(t) ~ {Dt)~ d / 2 . Hence it is 
related to the spectral properties of the transition rate 
matrix. 

In the d=l case, it is well known [20] that D can dis- 
play an abrupt percolation-like transition from diffusive 
(D > 0) to sub-diffusive (D = 0) behavior, as the sparsity 
parameter drops below the critical value s cr — I. Simi- 
lar anomalies are found for fractal structures with d < 2, 
also known as "random walk on percolating clusters" , 
see [2TH25] . A question arises whether such a transition 
might happen in higher dimensions. 

In [TUJ the spectral properties in the d=2 case were 
investigated: on the basis of the renormalization group 
(RG) procedure it was deduced that V(t) decays in a 
logarithmic way, indicating anomalous (sub) diffusion. In 
the present work we shall introduce a different approach 
that implies, contrary to the simple RG treatment, that 
in spite of the sparsity, the long time dynamics is in fact 
diffusive rather than sub-diffusive. 

Resistor network picture.- One can regard the p n in 
Eq. ([I]) as the charge in site n; each site is assumed to 
have unit capacitance; hence p n —Pm is the potential dif- 
ference; and w nm {p m —p n ) is the current from m to n. 
Accordingly Eq. ([TJ) can be regarded as the Kirchhoff 
equation of the circuit. While calculating D it is illu- 
minating to exploit the implied formal analogy with a 
resistor network calculation [T2] HH HE] US] • Namely, re- 
garding w nm as connectors, it follows that D is formally 
like conductivity. It follows that D[w] is in general a 
semi-linear function: 



D[Xw] = XD[w] 
D[w a + w b ] > D[w a ] +D[w b ] 



(8) 
(9) 



If the rates are induced by a driving source, the above su- 
per additivity implies that the analysis should go beyond 
the familiar framework of linear-response theory |27j . 

In this work we obtain an improved estimate for D 
that we call effective range hopping (ERH). Using this 
approach we show that in the d s becomes 

small, the functional D[w] exhibits a smooth crossover 
from "linear" behavior to "semi-linear" VRH-type depen- 
dence. Our approach is inspired by the resistor network 
picture of [7J IT2HT0I 127] . and leads in the appropriate 
limit to the well known Mott's variable range hopping 
(VRH) estimate for D. 

Outline.- We first describe some known results, and 
some additional numerical results, for the spectral prop- 
erties of d=\ and d=2 networks, and for the dependence 



of D on the sparsity. Then we show that an ERH pro- 
cedure is useful in describing the crossover from the lin- 
ear regime (no sparsity) to the semi-linear regime. In 
the latter regime a "resistor network" approach is es- 
sential, and the percolation threshold manifests itself in 
the calculation. Finally we demonstrate that the same 
ERH procedure can be applied in the case of a quasi- 
one-dimensional network that is described by a sparse 
banded random matrix. The latter is of relevance to pre- 
vious studies of energy absorption by a weakly chaotic 
system [27]. We conclude with a discussion and a short 
summary. 



II. THE RANDOM SITE HOPPING MODEL 

Consider a network that consists of sites that are dis- 
tributed in space, locations x n . With each bond nm we 
associate an activation energy e nm > 0, and assume 



w e 



-\x n -x m \/i 



Accordingly we have the identification 
B{r) = e~ r/s 



(10) 



(11) 



We note that in the traditional formulation of the Mott 
problem the "activation energies" are not due to some 
"barriers" , but are determined by the on-site binding en- 
ergies, namely e nm = \s n — £ m \/T, where T is the temper- 
ature. In this paper we treat the e nm as an uncorrelated 
random variable. 

The density of sites relative to some initial site is char- 
acterized by a joint distribution function 

p(r, e)drde = f(e)de, n d = 2, 2tt, 4tt (12) 

r o 

We distinguish between the "Mott hopping model" and 
the "degenerate hopping model" . Namely, 

/(e) = I Mott hopping model (13) 

/(e) = <5(e) Degenerate hopping model (f4) 

The normalization of /(e) as defined above fixes the value 
of the constant r d , which we regard as the "unit cell". In 
the numerics we set the units of distance such that ro = I . 

In the traditional formulation of the Mott problem it 
is assumed that mean level spacing within £ d is A^, such 
that the number of accessible sites is (cfe/Aj) (d?r/t t d ). 
By the convention of Eq. (I2| this implies that the unit 



cell dimension is temperature dependent 



~0 



T 



e 1 



[for Mott model] 



(15) 



We re-emphasize that the number of sites per unit vol- 
ume in the Mott problem is infinite, but effectively only 
~ T/A^r sites are accessible within £ d per attempted tran- 
sition. It is convenient to characterize a random site 
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model by a "sparsity" parameter that is defined as in 
Eq.Q. Accordingly 



ro 



T 

a; 



l/d 



[for Mott model] (16) 



We refer to a network as "sparse" if s <C 1. 

The lattice model with near-neighbor (n.n.) transi- 
tions is one of the most popular models in statistical 
mechanics: in particular the random walk problem on 
a lattice is a standard textbook example. If the rates 
are generated from a log-wide distribution, it can be re- 
garded as a variant of the random site hopping model. 
For details see Appendix |Aj In particular we note that 
the d=l version is formally equivalent: it does not matter 
whether the distribution of w is due to random distances 
r, or due to random activation energies e. 

Finally we note that a quasi-one-dimensional version of 
the random site model arises in the study of energy ab- 
sorption as explained in Appendix|Bj and later addressed 
m Section KD 



Regarded as a transport coefficient D relates the prob- 
ability current to the density gradient. This is known as 
Fick's law. From the discussion in the Introduction it fol- 
lows that D is like the conductivity of a resistor network, 
which relates the electrical current to the voltage differ- 
ence. Some further details on the practical calculation 
of the conductivity are presented in Appendix [EJ On the 
basis of this analogy it should be clear that D[w] is in 
general a semi-linear function of the rates, see Eq. @. 



IV. EXACT AND NUMERICAL RESULTS FOR 
THE d=l LATTICE MODEL 

In the case of a d=l lattice model with n.n. transitions 
it is natural to use the notation w n = ui nj „_i. Pointing 
out the analogy with adding connectors in series the ex- 
pression for D is 



D = 



= [«>1] 



1 



■ w 



(21) 



III. THE CHARACTERIZATION OF 
TRANSPORT 

The long time dynamics that takes place on the net- 
work is characterized by the spreading S(t), and by the 
survival probability V(t). If the system is diffusive, these 
functions have the following functional form: 



S(t) 
V(t) 



(r 2 ) t 



(2d)Dt 



(4-rrDt) 



d/2 



(17) 
(18) 



See Appendix [C] for details. The diffusion coefficient D 
appears here in consistency with its definition in Eq.([7|. 
We note that in the case of sub-diffusion 



S(t) cx t c 



[a < 1] 



(19) 



which implies by Eq.Q that D = 0. 

The spectrum of the matrix w consists of the trivial 
eigenvalue Ao = that is associated with a uniform distri- 
bution, and a set of negative numbers — A& that describe 
the decaying modes. The spectral function N{\) counts 
the number of eigenvalues up to the value A. We nor- 
malize it per site such that N{oo) = 1. The associated 
density of eigenvalues g(X) is related to Vit) by a Laplace 
transform. See Appendix [C] for details. It follows that in 
the case of a diffusive system 



g(X)d\ 



r 


'A" 




D 



d/2 



(20) 



In Appendix[D]we clarify that this expression agrees with 
Debye law. Accordingly the calculation of D parallels the 
calculation of the speed of sound c in Debye model. 



The calculation that leads to the last equality has been 
done with the distribution of Eq. ( |A2[ ), where s = £,/r$. 
Note that we have here a serial addition of resistors R = 
J2 n wnere Rn = l/w„. For s < 1 the distribution of 
each R n is dominated by the large values, hence R = oo. 
On the other extreme for s > 1 the distribution of the 
R n has finite first and second moments, and accordingly 
the result for R becomes self-averaging, as implied by the 
central limit theorem. This means the D is "well defined" 
only for s > 2. For 1 < s < 2 the result for the average 
R is finite but not self-averaging. 

The dependence of D on s is illustrated in FigJTJi. In 
the sub diffusive regime (s < 1), where the result for the 
diffusion coefficient is D = 0, the dynamics becomes sub- 
diffusive. The explicit results for the survival probability 
and for the spreading are known |20j : 



S(t) 
V(t) 



t 2s/(l + s) 
t -s/{l+s) 



and the associated spectral function is: 
Af(\) - A s/(1+s) 



(22) 
(23) 



(24) 



The numerical demonstration of the latter expectation 
is displayed in Fig.[2] (left upper panel). We clearly see 
that for s < 1 the asymptotic slope corresponds to sub- 
diffusion, while for s > 1 it corresponds to diffusion. 



V. NUMERICAL RESULTS FOR THE d=2 
RANDOM SITE MODEL 

Results for the spectral counting function of the degen- 
erate d=2 random site model are presented in Fig[2] (right 
upper panel). We also display there (in the lower panel) 
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the participation number (PN) for each eigenstate. The 
PN of an eigenstate that corresponds to an eigenvalue 
is conventionally defined as follows: 



PN = 



-, -l 



ElHA fe )l 4 



(25) 



As expected from the study of localization in a disor- 
dered elastic medium [28] , the PN becomes larger in the 
limit A — ► 0, without apparent indication for a mobility 
threshold. 

Assuming localized modes that are conceived via 
dimerization of neighboring sites, Af(\) should equal 
the probability exp[— V(r)/ro] not to have any neigh- 
boring site within the volume V(r) of the sphere 
2w exp(— r/£) > A. The RG analysis of [10 refines this 
naive expectation, adding a factor of 2 in the exponent, 
leading to 



Af(\) = exp 



n 



2d V \2w 



(26) 



where s = £/Vo. This expectation is represented in Fig[2] 
(right upper panel) by solid lines. We see that it fails 
to capture the small A regime, where the distribution 
corresponds to diffusive behavior. 

Extracting D via fitting to Eq.(20l we get Fig(T|3. We 
see that in the d=2 model there is no abrupt crossover 
to sub- diffusion. We therefore would like to find a way 
to calculate D, and hence to have the way to determine 
the small A asymptotics. 

Note added. - One should conclude that the RG of ref- 
erence [TU] applies only for the analysis of the high fre- 
quency response, while our interest is focused in the low 
frequency (dc) analysis. The crossover between the two 
regimes is implied. For more details in this direction see 
a follow-up work |29) that confirms our physical picture 
and demonstrates numerically the implied crossover. 



VI. LINEAR AND ERH ESTIMATES FOR THE 
DIFFUSION COEFFICIENT 



This expression is strictly linear. It describes correctly 
the average transient spreading. In the absence of disor- 
der we can trust it for arbitrary long time. But if we have 
a disordered or sparse network, the possibility for trans- 
port is related to the theory of percolation [71 SSI EH- We 
are therefore motivated to introduce an approximation 
scheme that takes the percolation aspect into account. 
We shall refer to this scheme as "effective range hop- 
ping" (ERH) because it is a variation on the well known 
VRH procedure. 

Inspired by [71[T21[T3] we look for the threshold w c that 
is required for percolation. In the ERH scheme we sug- 
gest using the following equation for its determination: 



|| p(r,e)drde 



(29) 



w(r 1 e)>w c 



Here n c is the effective coordination number that is re- 
quired for getting a connected sequences of transitions. 
For a d=2 square lattice model it is reasonable to set 
n c = 2, reflecting the idea of forming a simple chain of 
transitions. Rephrased differently the requirement is to 
have an average of 50% connecting bonds per site. For 
a d=2 random site model one should be familiar with 
the problem of percolation in a system that consists of 
randomly distributed discs. The effective coordination 
number that is required for getting percolation in such a 
model is n c — 4.5, as found in [30j . and further discussed 
in Section IV. A. 1 of [2]. 

The second step in the ERH scheme is to form an ef- 
fective network whose sparse elements are suppressed to 
the threshold value. Then it is possible to use the linear 



formula Eq. ( 28 1 . Hence we get 



^ERH — 



JJmin{w(r, e),w c } r 2 p(r, e) dedr 



(30) 



This expression, as required, is semi-linear rather than 
linear. It looks like the linear estimate of Eq. (28 1, but 
it involves a network with w nm that are equal or smaller 
to the original values. The "suppressed" connectors are 
those that are too sparse to form percolating trajectories. 



The standard way to calculate diffusion in a d=l ran- 
dom walk problem is to inspect the transient growth of 
the variance Var(n) = 2Dt. In the stochastic context, 
if we start at site n we have Var(n) = ^2 n , p n '{ n ' — n ) 2 , 
with p n i = w n i n t, hence 



D„ 



n) 2 w r . 



(27) 



The generalization to more than one dimension is 
straightforward. Averaging the transient expression over 
the starting point we get the result 



jTj JJ U '( r ' C ) ^ P( r ' 6 ) dedr 



(28) 



VII. VARIABLE RANGE HOPPING (VRH) 
ESTIMATE 



The ERH is similar to the generalized VRH procedure 
that we have used in previous publications [HI [19] . The 
traditional VRH is based on the idea of associating an 
energy cost e(r) to a jump that has range r. Namely, 



e(r) 



n 



A 



(31) 



corresponding to the average level spacing of the sites 
within a range r. In our notations e(r) = e(r)/T. For the 



5 



general network models that we consider here, the rela- 
tion between e and r is determined through the equation 



o Jo 



p(r',e')dr'de' = 



(32) 



where n* is of order unity. In fact we shall deduce later, 
in Section [Xj that for consistency with the ERH estimate 
this value should be n* = n c /d. With the substitution of 



Eq.(12| the trade-off equation can be written as 



F(e) 



(33) 



where F(e) is the cumulative distribution function that 
corresponds to the density /(e). In the Mott problem 
F(e) = e, and Eq.(31 1 is recovered. In words Eq.(32 1 asks 
what is the e window that is required in order to guaran- 
tee that the particle will be able to find with probability 
of order unity an accessible site within a range r. Larger 
jumps allow smaller cost. Then we estimate D as follows: 



w x 



(34) 



where r* is the optimal range that maximizes w(r,e(r)), 
with associated energy cost e* = e(r*), and effective tran- 
sition rate w* — w(r*,e*). See Figj3]for illustration. 

The VRH estimate, unlike the ERH, does not interpo- 
late with the linear regime. It can be used to estimate 
D only if the system is very sparse (s <C 1). It can be 
regarded as an asymptotic evaluation of the ERH inte- 
gral: it assumes that the hopping is dominated by the 
vicinity of the optimal point (r* , e* ) . Accordingly, VRH- 
to-ERH consistency requires the identification w* = w c . 
However, using known results from percolation theory, 
one possibly can further refine the determination of the 
optimal value w*. Namely, a somewhat smaller value 
than the threshold value w c might allow a better connec- 
tivity. As s becomes very small, the effective range Sw in 
the ERH integral, which contains the dominant contribu- 
tion, becomes very small compared with w c — w*, and one 
should be worried about the implied (sub-dominant) cor- 
rection. This speculative crossover is beyond the scope 
of the present study, and possibly very hard to detect 
numerically. A useful analogy here is with the crossover 
from "mean-field" to "critical" behavior in the theory of 
phase transition, as implied by the Ginzburg criterion. 



VIII. ERH CALCULATION FOR THE d=2 
LATTICE MODEL 

The d=2 lattice model, as defined in Appendix [A"} is the 
simplest and most common example for studies of perco- 
lation and percolation-related problems. We substitute 



ordination number Cj, 
the median value of the n.n 



into Eq. ( 29 1 the effective density Eq. ( A3 ) with the co- 
4, and deduce that w c is merely 
transition rates. The ERH 



calculation using Eq.(30l with Eq.(A3) requires a simple 



f(e)de integration, which can be rewritten as f(w)dw in 
tegral. This integral is the sum of w > w c and w < w, 
contributions, namely 

1 1 

2 



2 Wc 



wf(w)dw 



(35) 



Note that the first term in the square brackets originates 
from the w > w c contribution. Note also that the result 



fn for a delta distribution, i.e. in the absence 



is D = 
of disorder. 



IX. ERH CALCULATION FOR THE 
DEGENERATE HOPPING MODEL 

We now turn to the calculation of the ERH estimate 
for the degenerate hopping model. The ERH threshold 
can be written as w c = wq exp(— r c /£), where r c is deter- 



mined through Eq.(29), which takes the form 
Qcir^dr 



leading to 



wq exp 



a. 



l/d 



ro 



(36) 



(37) 



(38) 



The calculation of the ERH integral of Eq.( 30 ) is detailed 
in Appendix [F] We note that the linear approximation 
of Eq.(28) is formally obtained by setting r c = 0, leading 
to 

Au- = { -^^s^ WQ rl (39) 

Then it is possible to write the result of the ERH integral 

as 

£> EHH = EXP d+2 fy\ C- 1 /^ Aiaea, (40) 



where s c = £/r c , and 



EXP, (a;) 



i i 



k=0 



k\ 



(41) 



The linear result is formally obtained by setting n c = 
or in the d — > oo limit. In the other extreme of s <C 1 we 
get a VRH-like dependence 



D 



-l/s 



for s < f 



(42) 



Numerical verification.- To obtain an ERH estimate 
we have to fix the parameter n c in Eq.( 36 ) . One approach 
is to regard it as a free fitting parameter. But it is of 
course better not to use any fitting parameters. Fortu- 
nately we know from [30l|3T] that n c — 4.5 is the average 



number of bonds required to get percolation. The veri- 
fication of the ERH estimate for the random site model 
with this value is demonstrated in FigJIJb). 
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X. ERH CALCULATION FOR THE MOTT 
HOPPING MODEL 

We turn to calculating the ERH estimate for the non- 
degenerate Mott hopping model, and contrast it with the 
linear approximation, and with the traditional VRH esti- 



mate. The ERH threshold is determined through Eq.( 29 1, 
leading to 



w exp(-e c ) 



d n c 



i/(d+i) 



(43) 
(44) 



In the VRH procedure the optimal hopping range is 
found by maximizing w(r, e) along the trade-off line of 
Eq. (33 1, as illustrated in Fig(3j leading to 



i/(d+i) 



'0 



and the associated rate is 



provided n* = n c /d 



(45) 



(46) 



This identification is nexessary if we want the VRH to 
describe correctly the asymptotic dependence of D on s. 



The calculation of the ERH integral of Eq. ( 30 ) is de- 
tailed in Appendix [F] Thanks to our conventions the lin- 
ear result is the same as Eq.(39), and the the final result 
can be written as follows: 



= EXP 



d+3 



(e c ) 



Aine 



(47) 



where EXP(x) is the polynomial defined in Eq.(41). The 
linear result is formally obtained by setting e c = or in 
the d — > oo limit. 

We see that the VRH estimate can be regarded as an 
asymptotic approximation that holds for sCl. Using 
Eq. (151) and Eq. (16) we deduce from Eq. (39) and from 
Eq.(47| that 



while for s <C 1 



A 



A, 



2/(<2+l) 



T 



exp 



i/(d+iy 



(48) 



(49) 



where Tn is a constant. 



XL ERH CALCULATION FOR THE BANDED 
QUASI-ONE-DIMENSIONAL MODEL 

We can apply the ERH calculation also to the case of 
the quasi-one-dimensional model that we have studied in 
the past [HI Qjj] . This model is motivated by studies of 
energy absorption [27] ■ For details see Appendix [b| The 
network is defined by a banded matrix w. For simplicity 



we assume that the sites are equally spaced, and that the 
reason for the "sparsity" is the log-wide distribution of 
the in-band elements. 

The ERH threshold w c is deduced from Eq.fl29"j). For 
a general B(r) and /(e) one can integrated over de, and 
then it takes the form 



VLdr^dr 



F lO; 



w 



B(r) 



= n c (50) 



where F(e) is the cumulative distribution function that 
corresponds to the density /(e). Here we are considering 
a d=l network. However, we are dealing with a banded 
matrix which in some sense is like adding an extra (but 
bounded) dimension to the lattice. 

Specifically we assume that B(r) = 1 within the band, 
and zero for \r\ > b. The non zero elements have a log- 
box distribution, namely, e is distributed uniformly over 
a range [0, a]. To have large a means "sparsity". One 
should notice that this sparsity is less traumatic than 
having s <C 1 in the d—l lattice model that we have con- 
sidered in Section HVl This is because the distribution is 
bounded from below by a finite non-zero values. Accord- 
ingly we cannot have sub-diffusion here. 

We now turn to estimate D using the the ERH pro- 
cedure. It should be clear that the success here is not 
guaranteed for reasons that we further discuss in the last 
paragraph of this section. From Eq. (501 it follows that 
w c = wo exp(— e c ), where e c is the solution of 



26F(e c ) 



(51) 



For the assumed e distribution the solution of this equa- 
tion is trivial 



(52) 



While doing the ERH integral of Eq. ( 30 ) note that the 



integral dr should be replaced by a sum. It is convenient 
to define 



b = I> 2 = = \b(b + i)(2b + l) (53) 

r— 1 



Then the ERH estimate takes the form 
1 



D 



ERH 



1 



2b 



bw (54) 



The linear estimate of Eq. ( 28 1 is formally obtained by 



setting n c = 0, and in the absence of disorder it obviously 
reduced to D = bwo- We define 



9s = D/D 



linear 



(55) 



Numerical results are presented in Fig[4] and they agree 
with the ERH estimate. 

At this point one wonders whether D can be extracted 
from the spectral analysis, i.e. via fitting to Eq. (20 1. 
In Fig.[4]: we plot the " D" that is extracted from the 
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spectral analysis versus the D that has been found via 
the resistor network calculation. We observe that the 
obtained values are much smaller. Our interpretation for 
that is as follows: the density of eigenvalues is related 
to the survival probability V(t) via a Laplace transform; 
For a quasi-one-dimensional system there is a short time 
d=2 like relatively fast transient; Consequently the d—1 
decay holds only asymptotically with a smaller prefactor. 
Accordingly we do not know whether there is a wise way 
to deduce D from the spectral analysis in the case of a 
qausi-one-dimensional network. 

Concluding this section we would like to warn the 
reader that the use of the percolation picture in d=l is 
somewhat problematic: strictly speaking there is no per- 
colation transition. Obviously for 6=1 we are back with 
the d=l lattice model for which there is sub-diffusion if 
s < s cr with s cr = 1. However, if b is reasonably large, it 
is not feasible to encounter such an anomaly in practice. 
Even if the distribution is not bounded from below, the 
redundancy due to b > 1 would lower the effective value 
of s cr . Furthermore: in the Fermi-golden-rule picture 
(see next section) the occurrence of "weak links" along 
the band are practically not possible because the matrix 
elements V nm are not uncorrelated random variables. We 
can refer to this as the rigidity. This rigidity is implied 
by semi-classical considerations. 

XII. SEMILINEAR RESPONSE PERSPECTIVE 

Considering models of energy absorption, see (Ap- 
pendix IbJ, it is assumed that the transition rate w nm , 
between unperturbed energy levels m and n, is deter- 
mined by a driving source that has spectral content S(lo). 
The Fermi golden rule can be written as 

w nm = §(E n - E m ) \V nm \ 2 (56) 

where V nm is the perturbation matrix in the Hamilto- 
nian. Accordingly we can write instead of D = D[w] an 
implied relation D = D[S(lj)]. This relation is in gen- 
eral semi-linear. This means that only the first property 
below, which corresponds to Eq. ^ is satisfied, not the 
second one. 

D[XS(oj)] = \D[S{w)] (57) 
D[S a (u) + S b {u)] = D[S a (u)] +D[S b {u)] (58) 

To have a semilinear rather than linear response may 
serve as an experimental signature for the applicability 
of resistor-network modeling of energy absorption. We 
note, however, that if the the driving were added "on 
top" of a bath, the response would become linear at small 
intensities. Namely, if one substituted 

S(w) total = 5 bath (u;) + S(oj) (59) 

it would be possible to linearize D with respect to the 
S(uj) of the driving source. 



The statement that VRH is a "semilinear response" 
theory rather than "linear response" theory is a source 
for non- constructive debates on terminology. The reason 
for the confusion about this point is related to the phys- 
ical context. Do we calculate "current vs bias" or do we 
calculate "diffusion vs driving" . The response is linear in 
the former sense, but semi-linear in the latter sense. 



XIII. DISCUSSION 

It should be clear that there are two major routes in 
developing a theory for D. Instead of deducing it from 
spectral properties as in |10j . one can try to find ways 
to evaluate it directly via a resistor network calculation 
[T2TfT5"] , leading in the standard Mott problem to the 
VRH estimate for D. 

In [TS1 [Tni HZ] this approach was extended to handle 
"sparse" banded matrices whose elements have log-wide 
distribution, leading to a generalized VRH estimate. In 
this work we have pursued the same direction and ob- 
tained an improved estimate for D, the ERH estimate. 
Using this approach we showed that in the d=2 case, as 
s becomes small, the functional D[w] exhibits a smooth 
crossover from "linear" behavior to "semi-linear" VRH- 
type dependence. 

Relation to other models. - Disregarding the "sparsity" 
issue, the model that we were considering is a close rel- 
ative of the Anderson localization problem. However it 
is not the same problem, and there are important differ- 
ences that we would like to highlight. For the purpose of 
this discussion it is useful to be reminded that the hop- 
ping problem that we have addressed is essentially the 
same as studying the spectrum of vibrations in a disor- 
dered elastic medium. Hence D 1 / 2 parallels the speed of 
sound c of the Debye model. See Appendix [D] 

Mott vs Anderson.- In the hopping model all the off 
diagonal elements are positive numbers, while the nega- 
tive diagonal elements compensate them. It follows that 
we cannot have "destructive interference" , and therefore 
we do not have genuine Anderson localization. Conse- 
quently in general we might have diffusion, even in d=l. 
In d=2 we have a percolation threshold, which is again 
not like Anderson localization. See the discussion of frac- 
tons in [23] . 

Debye vs Anderson.- In the standard Anderson model 
the eigenvalues form a band A € [— A c , A c ]. The states at 
the edge of the band are always localized. The states in 
the middle of the band might be de- localized if d > 2. 
The spectrum that characterizes the hopping model does 
not have the same properties. With regard to the lo- 
calization of vibrations in a disordered elastic medium 
|28j . it has been found that the spectrum is A € [0, A c ]. 
The ground state is always the A = uniform state. The 
localization length diverges in the limit A — > 0. Conse- 
quently the Debye density of states is not violated: the 
spectrum is asymptotically the same as that of a diffu- 
sive (non-disordered) lattice. It follows that the survival 



probability should be like that of a diffusive system, and 
therefore we also expect, and get, diffusive behavior for 
the transport: spreading that obeys a diffusion equation. 



XIV. SUMMARY 

This was originally motivated by the necessity to im- 
prove the resistor-network analysis of the diffusion in 
quasi-one-dimensional networks [18 , and additionally 
from the desire to relate it to the recent RG studies [TU] 
of the spectral properties of random site networks. The 
key issue that we wanted to address was the crossover 
from linear-like to semi-linear dependence of D on the 
rates. This crossover show up as the "sparsity" of the 
system is varied. 

It should be clear that unlike the RG based expecta- 
tion of [TU], our analysis indicates that there is no sub- 
diffusive behavior in d—2. Accordingly, the anomalous 
logit) spreading that is predicted in [10] should be re- 
garded as a transient: for very small value of the sparsity 
parameter this transient might have a very long duration, 
but eventually normal diffusion takes over. 

One can regard "sparsity" as an extreme type of disor- 
der: the rates are distributed over many orders of mag- 
nitudes. Still, unlike the d=l case, the implication of 
"sparsity" in d=2 is not as dramatic: there is no "phase 
transition" between two different results, but a smooth 
crossover. It is therefore clear that our statements are 
consistent with those of older works that relate to the 
diverging localization properties of the low frequency vi- 
brations in disordered elastic medium |28j . 

The effective range hopping (ERH) procedure that we 
tested in this paper is a refinement of well known studies 
of variable range hopping ITlTfTS] . We used the in- 
sight of [13l Ej that connects VRH with the theory of 
percolation. 

Disregarding possible inaccuracy in the determination 
of the optimal rate, the ERH calculation provides a lower 
bound for D. Accordingly, by obtaining a non-zero result 
it is rigorously implied that D is finite. The purpose of 
the numerics was to demonstrate that in practice the 
outcome of the ERH calculation provides a very good 
estimate of the actual result, interpolating very well the 
departure from linearity. 

It was important for us to clarify that a large class of 
networks can be treated on an equal footing. In particular 
we demonstrated that the application of the ERH esti- 
mate does not require any fitting parameters. We have 
verified that the same prescription can be applied both in 
the case of the d=2 lattice model, and in the case of the 
d=2 random-site model, provided one uses the appropri- 
ate percolation threshold that is known from percolation 
theory. 

For the traditional Mott hopping model and its de- 
generated version we obtained the refined expressions 
Eq. ( 47 1 and Eq. ( 40 ) respectively. In these expressions 



and the crossover to linear response as a function of the 
sparsity (s) is transparent. Note that in the degener- 
ate random site model the sparsity is merely a geomet- 
rical feature, while in the non-degenerate Mott model 
the sparsity depends on the temperature as implied by 
Eq.Jlo). 



We would like to re-emphasize that the original moti- 
vation for this work is was the study of energy absorption 
by driven mesoscopic systems. In this context the impli- 
cation of the semi-linear crossover is the breakdown of 
linear response theory. The latter issue has been exten- 
sively discussed in past publications [27 . 
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Appendix A: Lattice model with n.n. hopping 

For s -C 1 the d= 1 random site model is essentially 
equivalent to a lattice model with equally spaced sites, 
near neighbor transitions, and random e. From the iden- 
tification e = r/£ it follows that the distribution of the 
"activation energy" is 

/(e) = s exp(-se), s = £/r (Al) 

This implies that the the distribution of the rates is 



f(w)dw = [w < Wo] 



(A2) 



The density of sites to which a transition can occur is 



p(r,e) = c L S{r - r ) /(e) 



(A3) 



where = 2 is the coordination number. This corre- 



sponds to the d=l case of Eq.(12). 

The d=2 version of the lattice model has no strict re- 
lation to the d=2 random site model. A popular choice 
is to assume a box distribution for the activation energy 
within some interval < e < a. The density of sites to 
which a transition can occur is 27rr/(e) for large r, as im- 
plied by Eq.(12). But for small r the effective density is 
given by Eq.(A3) with the coordination number cj, = 4. 



the full dependence on the dimensionality (d) is explicit, 



Appendix B: The quasi-one-dimensional banded 
matrix model 



On equal footing we consider the quasi-one- 
dimensional banded lattice model. This model is 
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motivated by studies of energy absorption [57]. In 
this context the transition rates are determined by the 
Fermi-Golden- Rule (FGR). Hence we write: 



w e 



B (E n — E„ 



(Bl) 



Here n and m are unperturbed energy levels of the sys- 
tem, but we shall keep calling them "sites" in order to 
avoid duplicated terminology. The density of sites rela- 
tive to some initial site is characterized by the same joint 
distribution function as for the d=l network, 



P(r, e) 



2/(e) 



(B2) 



Here r = \E n — E m \ is the distance between the energy 
levels, which is formally analogous to r — \x n — x m \ in the 
random site hopping model. We use here units such that 
the mean level spacing is unity. In the later numerical 
analysis we assume equally spaced levels such that the 
distance is simply r = \n — m\. 

In the physical context the band profile B(r) is deter- 
mined by the semiclassical limit, while the distribution 
of the e values is implied by the intensity statistics of the 
matrix elements. This intensity statistics is known as 
Porter-Thomas in the strongly chaotic case, correspond- 
ing to the Gaussian ensembles, but it becomes log-wide 
for systems with "weak quantum chaos" [32], reflecting 
the sparsity that shows up in the limiting case of inte- 
grable system |19j . 

In the numerical analysis we have considered simple 
banded matrices, for which B(r) = 1 for r < 6, and zero 
otherwise. Accordingly 1+26 is the bandwidth. The 
elements within the band are log-box distributed: this 
means that e is distributed uniformly over a range [0, <r]. 
Note that log-box distribution is typical of glassy sys- 
tems, where the tunneling rate depends exponentially on 
the distance between the sites. 



Appendix C: Numerical extraction of D 

In a diffusive system the coarse grained spreading is 
described by the standard diffusion equation, with an 
evolving Gaussian distribution 



1 



: exp 



2S x (t) 



(CI) 



fL\ y/2nS x (t) 

where S x (t) = 2Dt. It follows from this expression that 

S(t) = (r 2 {t)) = (2d)Dt (C2) 

Starting with all the probability concentrated in one 
"unit cell" we get for the survival probability 



Hi) 



(C3) 



(AnDt) d/2 

The eigenvalues of the diffusion equation are 

A fe = Dq 2 k , k = index (C4) 



where the possible values of the momentum are de- 
termined by the periodic boundary conditions as q — 
(2ir/L)k. It follows that the cumulative number of eigen- 
states per site is 



A/-(A) = 







' A" 


d/2 




d 


p 





(C5) 



It is well known that the survival probability is related 
to the eigenvalues of w through the relation 



-At _ 



j(X)dX e 



-At 



(C6) 



For a diffusive system one can verify that the expressions 
above for ^(A) and V(t) are indeed related by a Laplace 
transform. More generally, it follows that D can be de- 
duced from the asymptotic behavior of g(X) in the A — > 
limit where the diffusive description is valid. In contrast 
to that for large A we expect g(X) to coincide with the 
distribution of the decay rates j n = J2 m w mn, reflecting 
localized modes. 



Appendix D: Relation to Debye model 

Consider a system of units masses that are connected 
by springs. Once can describe the system by a ma- 
trix w whose of-diagonal elements w nm are the spring 
constants. The eigen-frequencies are determined accord- 
ingly, namely, u>k = \fXk- Assuming that the low lying 
modes are like acoustic phonons with dispersion u) = c\q\, 
where c is the so called speed of sound, one deduces that 



Wk\, 



k = index 



(Dl) 



Consequently the associated counting function is as in 
the Debye model: 



r Q \ d i)d 
d 



\2ncJ 



(D2) 



Co mpar ing the above expressions with Eq. ( C4 1 and 
Eq. ( C5 ) it follows that the calculation of c 2 is formally 



the same as the calculation of D. 



Appendix E: The resistor network calculation 

The diffusion coefficient D is formally like the calcula- 
tion of the conductivity of the network. Therefore it can 
be determined via a numerical solution of a circuit equa- 
tion. It is convenient to use the language of electrical 
engineering to explain how the resistor network calcu- 
lation is carried out in practice. Accordingly we use in 
this appendix the notation G instead of w for the matrix 
that describes the resistor network, and a instead of D 
for its conductivity. We define a vector V = {V n }, where 
V n is the voltage at node n, analogous to p n . We also 
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define a vector J = {/„} of injected currents. The Kirch- 
hoff equation Eq.|l]) for a steady state can be written as 
GV = 0. 

If the nodes were connected to external "reservoirs" 
the KirchhofT equation would takes the form GV = I. 
The matrix G has an eigenvalue zero which is associated 
with a uniform voltage eigenvector. Therefore, it has a 
pseudo-inverse rather than an inverse, and consequently 
the Kirchhoff equation has a solution if and only if the 
net current is J2 n I n — 0. 

For the purpose of calculating the conductivity we add 
a source 1\ = — 1 and a drain 1% = 1. We select the loca- 
tion of the source (site #1) and the drain (site #2) away 
from the endpoints. From the solution of the Kirchhoff 
equation we deduce 

[{V 2 -V 1 )/L}- 1 



a[d=l] 



(El) 



where L is the distance between the contacts. 

With regard to the quasi-one-dimensional model, we 
take the distance between the contacts to be L' = N/2 
and look at the voltage drop along an inner segment of 
length L = L' — 2b, to avoid the transients at the contact 
points. 

To find the conductivity in the d=2 case we select con- 
tacts points that have distance L ~ (N/2) 1 / 2 , and use 
the formula 



a[d=2] 



[(V2-Vx)/ln(L/e)] 



(E2) 



where I ~ 1 is the shift of the measurement point from 
the contact point. Here the voltage drop is divided by 
\n(L/£) instead of L, reflecting the two-dimensional ge- 
ometry of the flow. 



This leads directly to Eq.(40| with Eq.(39). 



Turning to the non-degenerated Mott model we have to 
deal with a two dimensional integral drde that has, as in 
the previous case, two domains w > w c and < w < r c . 
The two domains are separated by the line e + (r/£) = e c . 
It is therefore natural to change variables: 



+ (r/0 
H + (r/0) 



(F3) 
(F4) 



hence 



2dr\ 



1 1> />'- ( £ » +£ !)' 



2dr$ 
2dri 



d+l 



x/2 
c d+3 



2dr l (d + 2 



(d+2){d + 3) 
C d+2 T(d + 3,e c ) 



2d{d + 2)(d + 3)r l 



rfr(d + 4,e c ) 



(F5) 



This leads directly to Eq.(47| with Eq.(39). 



Appendix F: Calculation of the ERH integral 



The calculation of the ERH integral for the random 
site model involved the incomplete T function [33j . 



r(^+l, x)= ( r e c- r dr = i\ EXP 



t (x) e~ x (Fl) 



We first consider the degenerate Mott model. We sub- 
stitute in Eq.pOb, the w(r,e) of Eq.Q, and the p(r,e) 
of Eq.Q with Eq.Q. Thanks to the 6(e) we are left 
just with a dr integration that is split into the domains 
< r < r c and r > r c . Namely, 

r d+l 



^ERH — 



2d 

Wpttd 

2d 
2d 



,-ro/£'_ 



-dr 



r d 

' p 

e- r/i — ^dr 



'o 



-d+2 



d + 2 



wpn d z d+2 1 

2d 

wpttd£, 
2d(d + 2)r d * 



r d + 2. 



d+2 



r d + 3, 



(F2) 
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FIG. 1: Spreading in the d=l lattice model (a), and in the d—2 degenerate random site model (b). Panel (a) is based on 
known exact results. Its dashed blue line is the power a of the spreading, showing a sub-diffusive regime for s < 1, and a 
diffusive regime for s > 1. Its solid red line is the diffusion coefficient D, which is zero in the sub-diffusive regime. Panel (b) 
displays numerical results that refer to a network that consists of TV = 2000 sites randomly scattered over a square with periodic 
boundary conditions. The vertical axis is the diffusion coefficient D in a logarithmic scale, while the horizontal axis is X — —1/s. 
The numerical red dots are based on a resistor network calculation (see AppjEj, while the stars are extracted from the spectral 
analysis (see Fig|2|. The dashed line is the linear estimate (corresponds to n c = 0), while the solid line is the ERH estimate 
with n c = 4.5. One observes that the ERH calculation describes very well the departure from the linear prediction. 




FIG. 2: The cumulative eigenvalue distributions TV(A) for the d=l (ID) and for the d—2 (2D) models of Fig.[T] and the 
respective PN of the eigenstates (lower panels). Several representative values of s are considered. The dots are determined via 
numerical diagonalization of TV x TV matrices, each representing a network that consists of TV = 1000 sites randomly scattered 
over a square with periodic boundary conditions. There is a striking difference between the d—1 and the d=2 cases. For d=l, 
the log-log slope of A/"(A), see dashed lines, is less than d/2 for sparse networks (s < 1), meaning that we have sub-diffusion. 
In the d=2 case the small-A log-log slope is always d/2, which corresponds to normal diffusion. The solid lines in the upper 
2D plot are according to the RG analysis of [10) . namely Eq.(26l. The horizontal dashed line in the lower panels indicates the 
special value PN= 2 that corresponds to dimer formation. 



13 




FIG. 3: Comparing the VRH with the ERH procedure. The 
solid blue line that corresponds to the ERH threshold w c en- 
closes an "area" that corresponds to n c . The VRH trade-off 
is represented by the dashed red line. The VRH optimum 
is represented by the thic k re d dot. The VRH-to-ERH con- 
sistency requirement Eq. ( 46 1 is to have the VRH optimum 
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FIG. 4: We consider a quasi d—1 network that consists of 
N = 1000 sites with periodic boundary conditions. The net- 
work is described by a sparse banded matrix. The bandwidth 
is b, and the log-width of the rate distribution is a. See text 
for details, (a) The numerical result for g 8 = D/D Uaeax im- 
aged as a function of a and b. The values of D are found via a 
numerical resistor network calculation, see App(E] (b) Plot 
of the subset of results that refer to the b = 10 matrix. The 
curve is the ERH prediction, (c) Scatter diagram that shows 
the correlation between the " D" that is extracted from the 
spectral analysis, and the D that has been found via the re- 
sistor network calculation. 



